Quantum Monte Carlo study of a resonant Bose-Fermi mixture 
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We study a resonant Bose-Fermi mixture at zero temperature using the fixed-node Diffusion Monte 
Cario method. We expiore the system from weak to strong boson-fermion interaction, for different 
concentrations of the bosons relative to the fermion component. We focus on the case where the 
boson density ub is smaller than the fermion density np, for which a first-order quantum phase 
transition is found from a state with condensed bosons immersed in a Fermi sea, to a Fermi-Fermi 
mixture of composite fermions and unpaired fermions. We obtain the equation of state and the 
phase diagram and we find that the region of phase separation shrinks to zero for vanishing ub- 

PACS numbers: 67.85.Pq, 03.75.Ss, 03.75.Hh 
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Let us consider a system of bosons and spinlcss 
fermions with a tunable short-range boson-fermion (BF) 
attraction. For weak attraction, at sufficiently low tem- 
perature the bosons condense, while the fermions fill a 
Fermi sphere, and the BF interaction can be treated with 
perturbative methods [3, 0] • For sufficiently strong at- 
traction, bosons and fermions pair into molecules. In 
particular, for a fermion density np larger than the 
boson density tib one expects all the bosons to pair 
with fermions. The boson condensate is then absent in 
such a regime and the system should be described as 
a weakly interacting Fermi-Fermi mixture, one compo- 
nent consisting of molecules, with density tlm — tib, and 
the other component of unpaired fermions, with density 
njj = np — Ub- 

How does the system evolve at zero temperature be- 
tween the two above physical regimes? Several scenarios 
could be imagined in principle: (i) a continuous quantum- 
phase transition, with the condensate fraction vanishing 
smoothly at a certain critical value of the BF coupling; 
(ii) a first-order quantum phase transition, with phase- 
separation between a condensed phase and a molecular 
phase without condensate; (iii) the collapse of the system 
in the intermediate coupling region, with no stable state 
connecting the two different regimes. 

The above question has been object of increasing atten- 
tion recently, especially in the field of ultracold trapped 
gases, where the interaction can be tuned by using Fano- 
Feshbach resonances 0]. In particular, for "broad" res- 
onances, a Bose-Fermi mixture can be accurately de- 
scribed by a minimal set of parameters: the scatter- 
ing lengths cibb and cibf describing, respectively, the 
boson-boson and boson-fermion interaction; the boson 
and fermion densities ub and tif, and their masses tub 
and 77i p (the short-range fermion- fermion interaction be- 
ing negligible, due to Pauli exclusion). 

Initial experiments [!, 0] with ultracold Bose-Fermi 
mixtures supported the collapse scenario, with the insta- 
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FIG. 1. Energy of a BF mixture at x — 0.175 and f = 1, with 
the contribution of the bare binding energy of the molecules 
subtracted for a B F > 0. Circles: JS FN-DMC results. The 
solid line corresponds to eq. (|3]) and the dashed line is a guide 
to the eyes. Squares: JMS FN-DMC results. The dotted line 
corresponds to eq. ((5} with M* and a from and C = 0. 
Inset: Energy without subtracting the bare binding energy. 



bility occurring already for moderate BF coupling. How- 
ever, only a limited region of the parameter space was 
explored (e.g., a boson number Nb considerably greater 
than the fermion number Np and nonrcsonant values of 
the scattering lengths). The tunability of the Bose-Fermi 
interaction was first demonstrated in a 40 K— 87 Rb mix- 
ture , and then successfully exploited to form Feshbach 
molecules Recently an isotopic 40 K— 41 K mixture, 

exhibiting a broad Fano-Feshbach resonance to tune the 
BF interaction, was successfully cooled down to quantum 
degeneracy . The creation of Feshbach molecules has 
been finally achieved very recently also with 23 Na— 6 Li 
[lU and 23 Na- 40 K mixtures. The two above op- 
posite regimes of a Bose-Fermi mixture have thus been 
explored already to some extent in experiments. The in- 
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tcrmcdiate region, instead, has remained inaccessible so 
far, essentially because of the large losses due to three- 
body recombination processes onto deep energy levels. 
Some control of these losses should however be achieved 
by working with small concentrations of bosons, the dom- 
inant recombination process being proportional to n 2 B np. 
On the theoretical side, first studies 14-161 were 



strength is conveniently described in terms of the dimen- 
sionless parameter g = (kpaBF)^ 1 ■ 

For the BB repulsion, we set Rbb = l-086a BB and 
take C, = kpasB = 1- Such a constant BB repulsion guar- 
antees the stability for all considered values of the boson 
concentration x = ub/uf and of the BF couplings across 
the resonance. We notice that our value of ( is twice the 
critical value for mechanical stability found for g = in 
[23I ]. Due to computational time constraints, an analysis 
of the dependence on £ of the results is beyond the scope 
of this Letter. This question is however definitively rele- 
vant for experiments with ultracold atoms, since typical 
values of clbb and kp correspond to values of £ smaller 
by at least one order of magnitude than the value consid- 
ered here. In addition, with £ = 1, we expect the specific 
choice of the repulsive potential to play a role for 2; > 0.2, 
based on previous studies of bosonic systems (25. 26|. 

Simulations are carried out in a cubic box of vol- 
ume L 3 = Np/riF with periodic boundary condi- 
tions. We use a trial wavefunction of the general form 
V>t(R) = $s(R) ( E > a(R). ®s is a positive function of 
the particle coordinates R = (ri, . . . , rjv F , IV, . . . , rjv B ) 
and is symmetric under exchange of identical parti- 
cles, while &a satisfies the fermionic antisymmetry 
condition and determines the nodal surface of Vt- 
The symmetric part is chosen of the Jastrow form 
calculation (see Ref. for details). To our knowledge, $s(R) = Hi,i> fBF{riv)Y[ vj , /bb^.'j'), where the un- 

primed (primed) coordinates refer to fermions (bosons) 
and two-body correlation functions of the interparti- 
cle distance are introduced. In order to describe the 
Bose-Fermi and the Fermi-Fermi mixture, we use two 
different choices for $,4. The first choice (JS) is a 
usual Slater determinant for the bare fermions <&^(R) = 

A ^Vfei (l)V'fe2 (2)-"^fcjv F C^f)) i where A indicates the an- 
tisymmetrizer operator and ipk a {i) indicates a plane- wave 
state, with k a = 2n(n ax x + n ay y + n az z)/L and |A; Q | < 
kp. The second choice (JMS) is the antisymmetrized 
product of a Slater determinant for the molecules and a 
Slater determinant for the unpaired fermions 



the theoretical side, first studies 
based on mean-field or perturbative approaches and were 
mainly focused on the mechanical stability issue. The 
continuous quantum-phase transition scenario was first 
put forward in [l7| . within a "two-channel model" for 
the BF coupling (actually relevant for "narrow" Fano- 
Feshbach resonances [ 3. 1 1810 and further explored, for a 
broad resonance, in [19ll20|. The alternative scenario of a 
first order transition, with a rather vast region of phase- 
separation occurring between the "molecular" reg ion and 
the condensed one, was instead proposed in [2II |22|. A 
recent variational calculation has indicated finally that 
a sufficiently strong boson-boson (BB) repulsion should 
prevent the collapse scenario (23l |. 

In this Letter, we apply the fixed-node Diffusion 
Monte-Carlo (FN-DMC) method to the study of a res- 
onant Bose-Fermi mixture. This numerical technique 
yields an upper bound for the ground-state energy of the 
gas, resulting from an ansatz for the nodal surface of the 
many-body wavefunction that is kept fixed during the 



no theoretical work has exploited so far this technique to 
characterize the behavior of Bose-Fermi mixtures. We 
consider a three-dimensional homogeneous Bose-Fermi 
gas described by the Hamiltonian 
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+ J2 V B F(ru>)+J2 V BB(n> r ), (1) 

where ... and ... label, respectively, the fermions 
and the bosons. We consider equal masses raj = rap = 
m and model the interspecies BF interaction using an 
attractive square- well (SW) potential: VBF(r) = —V BF 
for r < Rbf (Vbf > 0), anc ^ zer0 otherwise, while the 
intra-species BB interaction is modeled by a repulsive 
soft-sphere (SS) potential: Vbb{t) = Vbb f° r r < Rbb 
(Vg B > 0), and zero otherwise. 

In order to eliminate any dependence on the range 
of the BF interaction potential we take Rbf such that 
npR BF = 10~ 7 or equivalently kpRsF = 0.0181 in 
terms of the Fermi wave vector kp = (d^np) 1 / 3 . In 
this regime the only dependence on the BF interaction 
potential is given by the scattering length o,bf- As is 
well known, clbf changes sign and diverges when a bound 
state appears (with binding energy eb = —h 2 /ma BF for 
o-bf > 0). In the many-body system, the BF coupling 
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where the molecular orbitals are defined as <fK a (i,i') = 
fb(\*i — exp (iK a (ri + iy)/2), which consist of the 
relative motion orbitals ft times the molecular center- 
of-mass plane- waves with \K a \ < Km, being um = 
K\,j/Q-k 2 } while for the unpaired fermions \k a \ < ku, be- 
ing 7i(7 = kfj/Qit 2 . The functions as well as Jbb and 
JbF: are chosen to be the solutions of the appropriate 
two-body problems, modified at long distance to com- 
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ply with periodic boundary conditions. Wc notice that 
the wavcfunction ([2]) is not symmetric under exchange of 
bosons. It is the analog of the Nosanow-Jastrow wave- 
function 271 \2Mi which has been successfully used in 



quantum Monte-Carlo studies of the equation of state 
of solid 4 He dilll. 

In Fig. Q] we report the FN-DMC results for the to- 
tal energy of the mixture at a small boson concentration 
x = 0.175 as a function of the interaction parameter g 
in units of the energy per particle of the free Fermi gas 
E-pQ = 3h 2 k F /10m = 3ef/5, where Ep is the Fermi en- 
ergy. We perform calculations with Np = 57, Nb = 10 
for the JS nodal surface ($^) and with N F = 40, 
Nb = 7 with the JMS wavcfunction ($^ /s ), in order 
to have almost equal bosonic concentrations. For the JS 
(JMS) nodal surface finite-size effects are considerably 
reduced by using closed shells for the number of fermions 
(molecules and unpaired fermions) and using Fermi liq- 
uid theory. The energy difference between the finite and 
infinite system, is assumed to be the same as for the non- 
interacting case (see [3lj|). We use this correction also to 
assess the error-bars, on top of the statistical error. 

In the weak-coupling limit we recover the perturba- 
tive results of Refs.[l|, |2|, which can be further expanded 
in powers of x, leading to the energy functional E = 
N f E fg (£ B f + Ebb), where 



£ 



BF 



1 



20 
9irg 



1 

1 + — 

■ng 



9ng 2 



1 



while Ebb = ( 1 + VxC 3/2 ^ 12 %r- ) corresponds to 

the energy of a weakly interacting Bose gas [32j. 

More generally, the condensed phase can be described 
in terms of a polaronic picture, where bosons are dressed 
by fermions. These polarons are characterized by an ef- 
fective binding energy A and an interaction term F. Sim- 
ilarly to [23[ one can thus introduce the following pola- 
ronic equation of state (EOS) holding in the limit i<1 



E pol = N F E FG [1 - A(g)x + F(g, ()x 2 ] 



(4) 



where A(g) = —(^r 
matrix framework [lj 



i)x^-o is calculated within a T- 
Hb being the chemical poten- 



tial of the bosons, while F(g, = 1&( 1 + 0)- An 
analogous (x 2 ) interaction term has been considered in 
the context of polarized Fermi gases [33 - l36j . In order to 
precisely evaluate the interaction coefficient D of the po- 
laronic EOS, we vary the concentration of the bosons in 
the relevant regime of couplings < g < 1 and we fit the 
coefficient directly from the FN-DMC simulations. Re- 
sults for the polaronic branch are shown in the inset of 
Fig. [2] in Table U we report the fitted values for D. The 
agreement with the polaronic EOS is rather good even at 
large concentrations. Some discrepancies start to appear 
at large x for g = 0.75. 

In Fig. [T] we compare the FN-DMC results with the 
JMS wavefunction to the energy functional E mo \ = 
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FIG. 2. Energy as a function of x, for ( = 1. Main figure, 
from top to bottom: energy of the FF mixture at g — 0.6, 
g = 0.75, g = 1 and g — 5, with the bare binding energy of 
the molecules subtracted. Dashed lines: eq. © with C = 0, 
M* = 2m and a = 1.18. Solid lines: best fit using eq. ([5]) 
with M* and a from Q]. Inset, from top to bottom: energy 
of the polaronic phase at g — 0, g = 0.5 and g = 0.75. Solid 
lines: best fit using eq. Q. 

NpE-pQ{~lQg 2 x/i + Eff), where the first contribution 
comes from the bare binding energy of the molecules and 



Eff = 



M*{g) 



[1 + xC(g, C)]x 5 / 3 + (1 - x) 5 / 3 



,, x 5a(o) ._. 
+ x(l-x)—^, (5) 
3ng 

which is expected to hold for large values of g. Here, 
the first term corresponds to the kinetic energy of the 
molecules whose effective mass isgiven by M*(g), taken 
from the analytic treatment of [J for a single molecule 
in a Fermi sea, corrected by a term proportional to the 
coefficient C for finite values of x. This higher order 
2-8/3 contribution could also embody a p-wave interaction 
between the molecules, which is expected to be significant 
for ( = 1 . The second and third terms correspond instead 
to the kinetic energy of the unpaired fermions and to 
the interaction energy of the two fermionic components 
which, at the level of mean-field theory, is proportional 
to the ratio a = dadfo-BF of the atom-dimer to the BF 
scattering length; the g— dependence of this coefficient is 
taken from [4| and in the strong coupling limit correctly 
reduces to the value a = 1.18 obtained from the solution 
of the three-body problem [53, [38| . At the small value 
of x = 0.175 shown in Fig. [Q the FN-DMC results well 
compare with the EOS $5$ with C = 0. 

Analogously to the polaronic branch, we perform sim- 
ulations using the JMS wavefunction for g > 0.6 and 
different concentrations of the bosons. Results of the 
molecular FF mixture are shown in Fig. [2] For the three 
largest values of g we find that the EOS in eq. (O, includ- 
ing the correction to M* linear in x, accurately describes 



TABLE I. Best fitted values of the D(g,Q coefficient of the 
polaronic EOS ([4]) and of the C(g, £) coefficient of the molec- 
ular EOS © for C = 1. 



D 



C 



0.00 0.99(1) 0.60 12.39(10) 
0.25 1.33(5) 0.75 3.37(2) 
0.50 1.75(5) 1.00 1.54(2) 
0.75 1.95(5) 5.00 0.60(1) 
1.00 1.25(1) 

the dependence of the energy on the concentration of 
bosons up to x = 1. The corresponding best fitted values 
of the coefficient C are reported in Table Q] For g = 0.6, 
our results start showing some deviations from the func- 
tional form ([5|) in the regime of intermediate concentra- 
tions 0.5 < x < 0.9. For even smaller values of g a num- 
ber of effects worsen the agreement with the FN-DMC 
data: i) the molecular effective mass from Q diverges for 
g ~ 0.5, indicating that a molecular picture is not valid 
anymore, ii) beyond mean-field interaction terms for the 
FF mixture are probably relevant and iii) the composite 
nature of the molecules should start to play a major role. 

We pass now to discuss the bosonic condensate fraction 
no = No /Nb j defined as the fraction of bosons occupying 
the k = single-particle state. For the polaronic phase 
one can determine no by calculating the unbiased long- 
range tail of the bosonic one-body density matrix from 
FN-DMC and Variational Monte Carlo (VMC) simula- 
tions (see e.g. Ref. [Sll]). The results for x = 0.175 are 
reported in the inset of Fig. [3] and show a constant de- 
crease of no from the weakly interacting regime, where 
the only contribution to depletion comes from the BB re- 
pulsion, to the strongly interacting regime, where the BF 
interaction dominates. In this region, however, the FN- 
DMC and VMC results for no start to differ significantly 
preventing an accurate unbiased determination of this 
quantity. As a consequence, we cannot assess whether 
no decreases to zero by following the polaronic branch 
deeper on the molecular side or, as it is more reasonable, 
remains always finite. In the molecular phase, which for 
large values of g is well described by the FF mixture of 
cq. (O, one expects no = 0. 

Fig. [3] finally presents the phase diagram in the x-g 
plane, for small concentration x, where we are confident 
in the validity of the energy functionals flU, © describing 
respectively the supcrfluid (SF) polaronic and the normal 
(N) FF molecular phase. The dotted curve corresponds 
to the energy-crossing between the two phases. The two 
homogeneous phases are separated by a phase separation 
region, obtained by finding the global minimum of the en- 
ergy functional v P E po \(gp,x P ) + (1 - v P )E mo \(g M , %) 
with respect to the fractional volume vp of the pola- 
ronic phase, and the local couplings and concentrations 
gp,xp,gM,XM, at fixed global particle numbers. (We 
have checked that, at equilibrium, the resulting local cou- 
plings and concentrations lie in the respective regions of 
validity of the two energy functionals.) 
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FIG. 3. Phase diagram in the x-g plane, at £ = 1. The 
dotted curve indicates the energy-crossing between the pola- 
ronic (SF) and FF molecular (N) phases. It lies in the mid- 
dle of a phase separation (PS) region. The gray area above 
x ~ 0.2 indicates the non-universal region (for £ = 1). Inset: 
bosonic condensate fraction no for x = 0.175 in the polaronic 
phase. Circles: extrapolated results from FN-DMC and VMC 
simulations. Dashed line: weakly-interacting Bose gas result 
Uq = 1 — (8/3)-^/a;C 3 /67r 3 . For g > 0.75 the unbiased estimate 
of no is no longer reliable resulting in large error bars. 



Our FN-DMC calculations support the scenario of 
a first-order quantum-phase transition, with a narrow 
phase separation region intervening between the con- 
densed polaronic and FF molecular phase. The width of 
the phase separation region shrinks to zero in the limit 
x — > 0, where the transition line tends to the critical 
coupling for the polaron-molecule transition, previo usly 
studied in the context of polarized Fermi gases [4, 4(| 
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However, contrary to Fermi gases, the polaron-molecule 
transition is not masked by a large phase separation re- 
gion with finite width even in the limit x — > [42] ■ The 
physical reason for this is the different contribution of 
molecules to the EOS: in polarized Fermi gases molecules 
are composite bosons and feel a repulsion mediated by 
the dimer-dimer scattering length, here they are compos- 
ite fermions and feel the Pauli repulsion. The polaron- 
molecule transition is thus connected continuously to a 
quantum phase transition occurring at finite boson con- 
centration, a result which does not depend on the value of 
£, provided the polaronic phase is mechanically stable, as 
can be seen from a small x expansion of the energy func- 
tionals. We note finally that such a conclusion would still 
be valid even if the first-order quantum phase transition 
found here was replaced by a continuous phase transi- 
tion, a scenario which we cannot completely exclude as 
our results depend on the choice of the competing states 
and on their EOS. 
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